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We apply black-box methods, i.e. where the performance of the method does not depend upon 
initial guesses, to extract excited-state energies from Euclidean-time hadron correlation functions. 
In particular, we extend the widely used effective-mass method to incorporate multiple correlation 
functions and produce effective mass estimates for multiple excited states. In general, these excited- 
state effective masses will be determined by finding the roots of some polynomial. We demonstrate 
the method using sample lattice data to determine excited-state energies of the nucleon and compare 
" the results to other energy-level finding techniques. 
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I. INTRODUCTION 



Lattice quantum chromodynamics (LQCD) has been used successfully to compute many experimentally observable 
quantities from first-principles calculation of Euclidean-time hadron correlation functions, even occasionally predicting 
experimental results before they are measured. However, the successes of LQCD have mostly been restricted to 
^ computation of the physical properties of the lowest-energy states in each quantum number channel by focusing on 
I ' the large-time behavior of correlation functions where uncertainties due to excited-state contributions are exponentially 
suppressed. Given that signal-to-noise in correlation functions also falls exponentially at large times, success is often 
(-H ' dictated by available computational resources. 

Both in meson and baryon spectroscopy there are many experimentally observed excited states whose physical 
properties are poorly understood that could use theoretical input from LQCD to solidify their identification. Other 
excited-state quantities that could be computed on the lattice, such as form factors and coupling constants, would 
be useful to groups such as the Excited Baryon Analysis Center (EBAC) at Jefferson Lab, where dynamical reaction 
models have been developed to interpret experimentally observed properties of excited nucleons in terms of QCD [ij, 0] ■ 
, In certain cases, input from the lattice may be helpful in determining the composition of controversial states, which 
^sj ' may be interpreted as ordinary hadrons, tetra- or pentaquarks, hadronic molecules or unbound resonances. 

. Among the excited nucleon states, the nature of the Roper resonance, iV(1440) Pn, has been the subject of interest 

■ since its discovery in the 1960's. It is quite surprising that the rest energy of the first excited state of the nucleon 
Q^ I is less than the ground-state energy of nucleon's negative-parity partner, the A^(1535)S'ii [3], a phenomenon never 

■ observed in meson systems. There are several interpretations of the Roper state, for example, as the hybrid state 
^ ' that couples predominantly to QCD currents with some gluonic contribution Q or as a five-quark (meson-baryon) 

iTj ' ^^^^^ 

Early LQCD calculations using the quenched approximation 0, H, M, ITol-Jlll [13 1 found the computed spectrum 
inverted relative to experiment, with Pn heavier than the Sn. A recent study [7| suggested that qualitative agreement 
between experiment and LQCD in the quenched approximation could be restored provided other simulation effects due 
to finite volumes and unphysically heavy quarks were properly addressed. The study strongly suggests the nature of 
the Roper resonance changes dramatically as the quarks are made physically light in LQCD simulations, as in Fig. [1] 
Clearly, future LQCD calculations will require improved analysis techniques for extracting multiple excited-state 
energies, as well as variational wavef unctions, in the nucleon sector to test the validity of this claim. 

Apart from the vast amount of detail about excited states accessible to LQCD computations with advanced analysis 
methods, the statistical accuracy of ground-state quantities is also enhanced because correlation functions computed at 
shorter Euclidean times can be used where the signal-to-noise is greater. As current lattice simulations are performed 
with ever-greater resolution at short Euclidean times as lattice spacings are decreased toward the continuum limit. 
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TABLE L Summary of existing published Sn and Pn calculations. Due to space limitations, we adopt those abbreviations for 
fermion actions: Domain- Wall Fermions [H, [13, [H, [ll (DWF), Chirally Improved Dirac Operator 20, 21] (CIDO), Fat-Link 
Irrelevant Clover 22] (FLIC); and for the analysis methods: Variational Method [23.[23| (VM). Constrained Curve Fitting (25j] 
(CCF), Maximum Entropy Method (MEM), Simplified Black Box Method (SBBM). For those works which do 

not perform extrapolation, we use the lightest pion mass to represent their results. 
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FIG. 1: Summary of previous lattice calculations with extrapolation to the physical pion mass point or the lowest simulated 
pion point (labeled as "f"). 

simultaneous extraction of ground and excited-state quantities will be essential to extract full value from such large- 
scale (and expensive) computations. 

A variety of analysis techniques have been applied to extracting the excited-state spectrum from correlation func- 
tions. The most widely used method is a nonlinear least-squares (NLLS) fit to a model function, such as a sum over 
two or more exponentials. Operationally, even such a simple nonlinear fit can be fraught with difficulty, from estab- 
lishing the range of Euclidean times included in the dataset to stabilizing the convergence of minimization algorithms 
by careful choices of initial guesses or temporarily freezing selected fit parameters during the minimization process, 
all of which require intervention by a trained expert. 

At such times, the expert typically turns to black-box methods for guidance because they do not require intervention 
to determine initial guesses and fitting ranges: estimates of correlation functions go into the black box and estimates 
of hadron energies come out. The main detraction of black-box methods are that the produced estimates are expected 
to have larger uncertainties, making them sub-optimal relative to least squares methods [29| . Marrying the two 
approaches can effectively combine the best features of both, leading to a highly-automated analysis program producing 
optimal estimates of energies. Prior to this work, black-box methods were mainly useful for extracting ground-state 
and, perhaps, first excited-state energies [1, [2^. We believe the method described below will provide the needed 
black-box method for estimating as many energies from fixed set of correlations functions as are likely to be extracted 
from a NLLS fit. 

The structure of this paper is as follows: Sec. [Til provides theoretical formulations of the excited-state effective 
masses for single and multiple correlators; it also explores certain extensions to these techniques: linear prediction 
and periodic boundary conditions. In Sec. IIIIl we apply the methodology to some characteristic lattice correlators 
and compare the results with simple fitting and the variational method [2j, |2^ . Conclusions and future outlook are 
given in Sec. IIVI Numerous details and examples are included in the appendices. Preliminary details of this work 
were presented in Ref. (soj . 
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II. THEORETICAL BASIS 



The hadron spectrum can be calculated in LQCD using two-point hadronic correlation functions 

C{to,t)^{o\0{t) O\to)\0), (1) 

where the creation and annihilation operators and O transform irrcducibly under the symmetries of the lattice 
space group [Sll. [ssl. [33| . After taking momentum (and spin for baryons) projection and inserting a complete set 
of hadronic eigenstates of the Hamiltonian (ignoring the variety of boundary condition choices possible) , this becomes 

M 

C{p, tn) = ^ A™(p) exp [-{to + na)E^{p)] (2) 

m— 1 

n > 0, A^, EmeR, 0<Ei<E2<---< Em- 

where Am contains not only the overlap factor between the eigenstate and states created by the operators but also 
any kinetic factors that do not depend on Euclidean time separation tn = na = t — t^. 



A. Effective Masses 



In general, a two-point correlation function computed on iV = 2M time-slices tn will admit an exact algebraic 
solution having the form of Eq. ^ with M energies E^ and amplitudes Am- The problem to solve is the nonlinear 
system of equations y = V(a:) a 
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(3) 



for Xm — exp [—aEm {p}] and — Am (p) exp [—toEm {p}] where yn — C {p, tn). V(a:;) is known as 2M xM rectangular 
Vandermonde matrix. 

By inspection, it appears the problem is of polynomial degree 2M and thus by the Abel-Ruffini theorem [H, [sH] 
should not admit a general closed form solution in terms of radicals for M > 2. The M = 1 solution is simple to 
compute and is widely known in the lattice QCD literature as the effective mass solution. Note already that the 
simple effective mass problem is linear and has only one solution, suggesting that the polynomial degree is actually 
of order M. 

The M = 2 solution was explicitly constructed by one of the authors [28J and was independently constructed some 
time later by others 8]. It was noted [28j that the problem, when reduced, required only the solution of a quadratic 
equation and so it was conjectured that the general problem of size M could be reduced to a polynomial equation in 
one variable of degree AI. 

An efficient algorithm has been available for some time for solving square Vandermonde systems [sj by making them 
upper triangular. This approach works equally well for rectangular Vandermonde systems as in Eq. ([3]). Furthermore, 
this approach reveals why the solution for the energies Em can be found without solving for the amplitudes Am and 
why the problem is of polynomial degree M. 

As a first step toward extracting the energies E from our data, we transform the system so that V{x) is in upper 
triangular form [sj by pre-multiplying by the lower 2M x 2M bi-diagonal matrices: 

1 
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(4) 
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where the first —1 on the diagonal appears in the m + 1 row and cohimn. In the appendices, we demonstrate in detail 
how the general solutions for M = 2, 3 and 4 work. 

Finding a general approach for M > 4 would be a tough challenge. Although Abel's Impossibility Theorem proves 
there are no general solutions in radicals for polynomials higher than quartic order, there arc numerical methods for 
finding the roots of polynomials of any order. The general form for the polynomial follows from Eqs. (jA10[) . (|A22[) 
and (|A39|) in the appendices: 
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Prony [38| showed that problems in the form of Eq. ([2]) implied the following system of equations y = H(y) p 
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(6) 



where the AI x M matrix H(y) has the special structure of a Hankel matrix and the components Pm of p are the 
coefficients of a polynomial 



M M 
m— 1 rn— 1 



(7) 



The Prony- Yule- Walker method (or just Prony's method, for short) [33, l39|, |40| solves Eq. ^ to find the coefficients 
and then finds the M roots of the polynomial in Eq. ([7]). The amplitudes are determined by substituting the roots 
into Eq. ^ and solving it. Note again that using 2 A/ timeslices of correlation function data to determine M effective 
masses is a problem of polynomial order M . 

The general conditions under which the solutions of the Hankel and Vandermonde systems coincide is presented in 
Ref. [4l| . Here we provide a simple demonstration that both solutions are the same under the assumption that there 
are no complications like degeneracies in the energy spectrum of Eq. Assuming H(?/) is invertible, solving Eq. ^ 
gives 



p = H-V, P(a:) = l + p^x=l + (H-V) x 

,2 



(8) 



for the polynomial of Eq. ([7]) and where x^ ~ ( ) . Recall that the inverse can be written in terms of the 

adjoint, or matrix of cofactors, H^^ = C/ |H|, Cij — (—1)*+-' |H(i; where the notation H(i; j) means removing 
row i and column j. For Prony's method, we can rescale P{x) — > |H| P{x) and still find the roots Xm by solving 



|H| + (Cy)' x = 



(9) 



for X. 



Returning to the Vandermonde method, the determinant of Eq. (0) can be expanded in terms of its cofactors 
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(10) 



\n\ = i-lf' |H| + J2 C.+i,M+ia:^ - (-1)''" |H1 + ^(-1)*+*^ \nii + 1;M+ l)\x\ 
As usual, each cofactor can be expanded in terms of further cofactors where additional rows and columns are removed: 

M 



\n{i + 1; M + 1)1 = ;^(-l)^+^ \n{l,i + M + 1)| y,. 



(11) 
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By eliminating the first row and last row of in Eq. ^ we recover H = ?<(l;Ai" + 1) and for the cofactors 

(-1^+1 \n{l,t + l;j,M + 1)1 = (-!)■'" |H(*;j)| (12) 

so the desired identity is recovered 

M M 

\n\ ^ i-ir |H| + i-lf'J2T.(-^y^' |H(*;j)|y,x\ (13) 

i=i ]=i 

up to a possible overall minus sign for odd M , which is irrelevant for finding roots. There is a unique set of solutions to 
the Vandermonde and Hankel systems (under the assumption of noise-free correlation functions with non-degenerate 
energy levels), so other considerations should determine which is the better method to construct the polynomial. It 
is our experience that computing coefficients from Eq. ([5]) is preferred to solving Eq. ([6]) as statistical noise in the 
correlation functions can lead to nearly singular Hankel matrices which are difhcult to invert. 

In an earlier work [2^, one of the authors showed that Prony's method (also called linear prediction) could easily 
be extended to use more than 2M timeslices of a correlation function to extract only M masses by constructing an 
over-constrained system of equations analogous to Eq. It is not obvious how to construct and solve a similar 
over-constrained system in the Vandermonde case. Thus, Prony's method has a potential advantage that more time 
samples of the correlation function can be used to extract the same number of energy levels leading to reduced 
statistical fluctuations. 



B. Solutions with multiple correlation functions 

When constructing correlation functions in LQCD, care is taken to ensure that the correlation function transforms 
irreducibly under the symmetries of the lattice space group [3l|, [H, [H, [s^l • For the model function, as in Eq. , this 
implies that the amplitudes depend on the details of the specific correlation function but that the energies depend 
only on the irreducible representation. Since it is common in lattice QCD simulations to compute at least two distinct 
correlation functions for each irreducible representation, effective mass solutions which combine data from multiple 
correlations are also possible. 

Assume that there are K correlation functions available as in Eq. ^ that differ only in their amplitudes: 

M 

Ck{p, tn) = ^ Akrnip) cxp [-{to + na)E^{p)] (14) 

m— 1 

n>0, Akrr.,EmeR, < Ei < E2 < ■ ■ ■ < E M ■ 

Data from the same TV time slices will be used in the following from each correlation function to construct M 
effective masses. Under this assumption the condition that there will be equal number of data points as unknowns is 
KN = [K + 1)M . In Appendix \K\ are the three solutions that satisfy the condition for i^T = 1, up to quartic order. 
There are four more solutions (up to quartic order) for K > 1: {K,M,N) = (2,2,3), (2,4,6), (3,3,4) and (4,4,5) 
(demonstrated in Appendix [B| . Relaxing the assumption that the same number of time slices are used from each 
correlation function will allow for more possibilities up to quartic order. It is straightforward to generalize to these 
cases if desired. 

The general form of the polynomial equation can be inferred by studying the solved examples in Eqs. (|B6[) . (jB16[) . 
(IB23[) and (|B29p . Define K Hankel matrices H^^^'^'' for each of the correlation functions with the constraints 
SitLi -^^fc = M and = M -I- 1. Then the general form of the polynomial equation is 
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As previously discussed, each Hankel matrix is generally of full column rank and, if the correlation functions are 
linearly independent, then the columns of different Hankel matrices are also linearly independent. So, Eq. (llSp will 
only be satisfied for discrete values of xi corresponding to the roots of the polynomial. 
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C. Periodic Boundary Conditions 

In practical LQCD calculations, the temporal extent is finite so the choice of temporal boundary conditions affects 
hadronic correlation functions near the boundary. For simplicity, starting from Eq. (U), set tQ = and identify the 
points to = and tj^ = Na which can be done using modular arithmetic, i.e. t„ = (n mod A'^)a. For anti-periodic 
boundary conditions, the typical hadronic Euclidean time correlation function is described by the model function 

M 

C(p,t„)= J2 { Arn^ exp [- {n mod N)aE,^{p)] + {-1 fAl{p)exp[~{N^n mod N)aE:jp)]} (16) 

m— 1 

n>0, A^, Al, E,n, E*^eR, < E^ < E2 < ■ ■ ■ < Em, < E*, < e; < ■ ■ ■ < EIj. 

For periodic boundary conditions, set (—1)^ 1. For mesons, B — but more importantly A,n — A^ and ~ E*^, 
which is not true for baryons {B — 1). So, baryon correlation functions represent M states propagating to the left 
and M different states propagating to the right for a total of 2M states. 

Meson correlation functions represent the same M states propagating to the right and left. However, time-reversal 
symmetry requires C(P,t„) — C{P,t]y^n), up to noise terms, so that only half of the computed timeslices are truly 
independent. Thus, as was the case with baryons, information about only M states in any given quantum number 
channel can be extracted from a single correlation function computed on 2M timeslices in a finite box. As shown in 
Ref. [28| . this can be made explicit by writing the meson correlation function as 

M 

C(r„) = ^ A™ exp{-aNE^/2) cosh{anE^), = {n - N/2)a. (17) 

m— 1 

To write this result in the Vandermonde form of Eq. ([3]), define the variables 

am = Amexp{~aNEm/2), Xm, ^ cosh{aEm,) , y-n = ^— r I . I C(t„_2j-i)- (18) 

j=o \ ^ J 

When solving Eq. Q, the domain of the solutions Xm will be the real numbers or complex conjugate pairs since 
real- valued correlation functions are used as input. Complex- valued solutions are clearly unphysical and should be 
discarded as they are likely due to noise. Real solutions may also be unphysical if they cannot be used to extract a 
non-negative energy, and this will depend on the details of the model function and hence the boundary conditions. 
For example, for the basic model of Eq. only the solutions < Xm < 1 will yield non-negative energies. For 
mesons in periodic boxes, Xm — cosh{aEm) so only Xm > 1 will yield non- negative energies. Finally, for baryons in 
periodic boxes, the M states propagating to the right have Xm = exp(— aE"™) and the M states propagating to the 
left have XM+m = exp(ai?^) so all solutions Xm > are physical and Xm > 1 means the state is propagating to the 
left. 

For some lattice fermion actions, e.g. staggered [13,111,1131 or domain- wall fermions ^4^,(161, a variation of Eq. ^ 
is needed as a starting point to account for states which oscillate in time. For example, an appropriate model for 
staggered mesons on an infinite lattice is 

M 

C{p, i„) - i-^^Em{p)] + {-irA*„M exp [-naE^M]} (19) 

m— 1 

n>0, Arn, A:^, Em, E*^eR, < Ei < E^ < ■ ■ ■ < Em, < E* < E* < ■ ■ ■ < El,. 

There are two independently ordered sets of states, half of which oscillate as (—1)". When solving Eq. ([3]), such 
oscillating solutions will have physical solutions if — 1 < Xm < 0. Similarly, for staggered baryons with periodic 
boundary conditions, physical solutions with Xm £ —1 are certainly expected as oscillating states moving to the left. 

Generally speaking, model functions appropriate for the common lattice discretizations and choice of boundary 
conditions can be formulated and rewritten in the Vandermonde form of Eq. ([3]) . The physical interpretation of the 
solutions Xm depends on the details of the discretization and boundary conditions. In some cases, all real solutions 
may have a physical interpretation and thus cannot be immediately discarded without further statistical analysis. 
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FIG. 2: EfTective mass plots from the 7 smeared-point proton correlators used in this work. The horizontal axis shows time 
and the vertical axis shows the effective masses, both in lattice units. 

III. NUMERICAL RESULTS 

Although the ground-state effective mass has a long tradition of use in the lattice QCD community, not much work 
has examined excited states and fewer yet their effective masses. In this section, we will demonstrate the application 
of these effective-mass techniques to some typical lattice correlation functions. At the end, we compare the results 
with the ones from variational method. 

The data on which we demonstrate these methods is from a study done using the quenched approximation to QCD, 
i.e. where the effects of quantum fluctuations of quark-antiquark pairs in the vacuum are ignored, greatly reducing 
the computational cost but leading to an unknown, but hopefully small, systematic error. We generated an ensemble 
of 16'^ X 64 anisotropic lattices with Wilson gauge action and nonperturbative clover fermion action using Dirichlet 
boundary condition. The spatial lattice spacing is about 0.125 fm with anisotropy 3 (that is, temporal spacing 
a^^ ~ 6 GeV). The parameters used in the fermion action give a 720-MeV pion. Specifically, our data are proton 
correlators using 7 Gaussian smearing parameters, 0.5-6.5 in steps of 1.0, including both smeared-point and smeared- 
smeared source-sink combinations. Fig. [2] shows the single-state effective-mass plots for these smeared-point proton 
correlators. 



A. Excited-Effective Masses 

We now apply the excited-effective masses to our nucleon data. In Fig. [3] we show the results of applying the 
single-correlator {K — 1) excited-effective masses to the nucleon data for all M with analytic solutions in terms of 
radicals. Notice that as the number of states included increases, the amount of early-time contamination is decreased. 
As the formulae better account for the exact form of the correlator, more of the time range can be reasonably used 
to determine the states. However, as the number of roots increases, the occurrence of "bad" roots (those that are 
negative or imaginary) tends to increase as well. Since these have to be thrown out, this causes gaps in the extracted 
states where the results are unreliable. 

In Fig. d] we show the results of applying the multiple correlator excited-effective masses to the nucleon data for 
all combinations of K and M that have solutions in terms of radicals. Notice that as the number of correlators 
increases, the quality of the extracted masses improves. The gaps in the data where bad roots appear become much 
less noticeable. 

We demonstrate this approach using the smallest smearing parameter (0.5) smeared-point correlator. There are a 
few parameters in the linear prediction method which we can tunc: the number of desired states L, the number of 
time slices used to predict the later time point N, and the order of the polynomial M. In this work, we will show a 
selection of the better choices in these degrees of freedom. Figure [S] shows the effective mass plot for L = 2,3,4 from a 
single Gaussian smeared-point correlator with fixed parameters = 20 and M — 8. The excited states are consistent 
with each other as one increases the value of K. Since we have used a large value of N to form the polynomial, 
each point uses information extracted from 20 timeslices. Thus, one does not need a large plateau to determine the 
final mass. One also notes that since we only use a single correlator to extract multiple states, the multiple states 
will be correlated; that is, large errors on higher-excited states will make the ground state noisy as well. A future 
improvement would naturally be to extend this approach to multiple correlators. 
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FIG. 3: Higher-effective mass plots with (top-to-bottom) one, two, three and four masses. The colors indicate the (black) 
ground, (red) first-excited, (blue) second-excited and (green) third-excited states. 



B. Variational Method 

Variational method [l^, [131 is a powerful tool for extracting multi-excited states in lattice QCD. We construct an 
r X r spectrum correlation matrix, Cij (t) , where each element of the matrix is a correlator composed from different 
smeared sources or operators Oi and Oj. Then we consider the generalized eigenvalue problem 

C{t)^j = X{t,to)C{h)^, (20) 

where determines the range of validity of our extraction of the lowest r eigenstates. If is too large, the highest- 
lying states will have exponentially decreased too far to have good signal-to-noise ratio; if is too small, many states 
above the r we can determine will contaminate our extraction. Over some intermediate range in to, we should find 
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FIG. 4: Multiple-correlator higher-eflFective mass plots with (top-to-bottom) {K,M} equal to {2,2}, {3,3}, {4,2} and {4,4}. 
The colors indicate the (black) ground, (red) first-excited, (blue) second-excited and (green) third-excited states. 

consistent results. 

If the eigenvector for this system is \a), and a goes from 1 to r. Thus the correlation matrix can be approximated 

as 



-tE„ 



(21) 



n=l 



with eigenvalues 



A„(<,io) = e-(*-*'')^" 



(22) 
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FIG. 5: Effective-mass plots from the linear prediction black-box method at fixed parameters N = 20 and M 



by solving 

Cito)-'^^C{t)Citor^/^ = \{t, to)^p. (23) 

Further analysis on the principal correlators, A„(f, io)j reveals information on the energy levels, En- 

The results from the linear prediction approach in Sec. Ill A] are compared with the variational method (4x4 with 
smearing parameter ranging 0.5-3.5), as shown in Fig. [S] Here we shift time with respect to the linear prediction 
plot by 10 to have a better comparison with the plateau region from the variational method. For the ground state, 
the numbers are consistent with the result from the variational approach, including the size of the error bar. This 
is remarkable, given that the amount of input information is a factor of 16 less in the linear prediction approach. 
The first-excited state is consistent but has larger error bar, which is no surprise. As for the second-excited state, 
it seems to be consistent with the variational ones but definitely needs more statistics. The third-excited state is 
much larger than expected from the variational approach, which might be caused by contamination from even higher 
excited states. 
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FIG. 6: Comparison between the variational method (4 x 4) and linear prediction with 4 extracted states 



IV. CONCLUSIONS 



The determination of the physical properties of the excited-state hadrons is currently of great interest due to the 
construction of the 12 GeV upgrade at Jefferson Lab, where such properties will be measured experimentally. Lattice 
QCD methods have the potential to predict these properties, provided they can be extracted from the exponential time 
series, called correlation functions, computed in Monte Carlo simulations. In this work, we demonstrate a powerful, 
yet easy to use, black-box method for analyzing one or more correlation functions and extracting information about 
excited states. It can easily be adapted to various choices of boundary conditions and discretizations. While the 
method can be used by itself to estimate physical properties of excited hadrons, we anticipate that it will also be 
useful as a method for generating initial guesses for nonlinear least-squares minimizcrs. 
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APPENDIX A: GENERAL EFFECTIVE MASS SOLUTION FOR M = 2, 3, AND 4 

1. General solution for M = 2 

To reduce the M = 2 problem, we pre-multiply by two factors of the bi-diagonal matrices of Eq. Q to find the 
reduced equation = L2LiVa. By introducing the auxiliary quantities: 

a. = xiy^-i ~y^ (2 < i < 2M) (Al) 
/3j = X2aj-i - aj (3 < j < 2M) (A2) 



the reduced system becomes: 



yi = ai+ 02 (A3) 

a2 = a2 (a;i - X2) (A4) 

P3 = (A5) 

(34, ^ (A6) 
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The first half of the equations, Eqs. (jA3p and (|A4p . involve both the amplitudes Ai,A2 and the energies Ei,E2 but 
the second half involve only the energies. It will be true for any M, in general, that the last M equations can be solved 
first to find all the energies. Once all the energies are known, the first M equations form a square upper triangular 
system that can be solved efficiently by backward substitution to find the amplitudes. 

To see that Eqs. (jA3p - (|A6|) yield the known solution [28,], first substitute Eq. (|A2|) and eliminate X2 from Eqs. (|A5p - 
dM} to find 



a2a4 — = 0, or 



Q!2 as 
as "4 



= 



(A7) 



where we note that the l.h.s. is the determinant of a 2 x 2 Hankel matrix or perhaps the minor of a larger Hankel 
matrix. After substituting Eq. (|Aip . this gives the known quadratic equation 

{vl - ym) xj + {yiVi - 2/2^3) xi + (^3 - y2yi) = 0. (A8) 
Note that this can also be written 



yi y2 

y2 ys 



yi 2/2 

ys yA 



Xi 



2/2 2/3 
2/3 2/4 



= 



(A9) 



where the coefficients are not determinants of Hankel matrices but minors of a single Hankel matrix. So, it can be 
written even more compactly as 



2/1 2/2 

2/2 2/3 
2/3 2/4 



(AlO) 



where the left block is a Hankel matrix and the right block is a Vandermonde matrix. 



General solution for M 



Using the auxiliary quantities defined in Eqs. (|Aip and (|A2p and a third auxiliary quantity: 

7, = X3A_i -P, (4 < ^ < 2M) 
the reduced system of equations for i\/ = 3 is: 



2/1 
a2 

74 
75 
76 



ai + 02 + 03 

02 (xi - X2) + a3 {xi ~ x^) 

03 {X2 ~ 2^3) {Xl - Xz) 







Following the procedure of Sec. lA 11 substitute Eq. (jAlip into the last three equations to eliminate x^ 
(redundant) set of equations 

Psk - Pi = 0, PsPe - P4P5 = 0, /?4/35 -Pl^O 



(All) 

(A12) 
(A13) 
(A14) 
(A15) 
(A16) 
(A17) 

and find the 
(A18) 



or equivalently 



/33 


Pa 


= 0, 


Ps 


Pa 


= 0, 


Pa 


P5 


Pa 


P5 




P5 


P6 


P5 


P6 



0. 



Note the l.h.s. of these three equations are the same as the three coefficients of Eq. (|A8p after the 
2/i Pi+2 and thus are minors of a Hankel matrix. Next, substitute Eq. (|A2p into these equations to 
and find the equation 



a2a4aQ + la-^a/^a^ — a| — a\a^ — a.2o\ = or 



a2 a-i 
as a4 as 

Q!4 Q!5 Qfg 



0. 



(A19) 

substitution 
eliminate X2 



(A20) 
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Finally, substituting Eq. (|Aip will produce a cubic equation in xi: 



Axl + Bxl+Cxi+ D = 

^ = 2/3 - 22/22/32/4 

B ^ 



(A21) 



yivl 



2/12/5 



2/32/4 + 2/22/4 + 2/22/32/5 
C = 2/32/4 - 2/32/5 - 2/22/42/5 -1 

2/22/5 



2/i2/6 



2/12/32/5 
2/12/42/5 

2/22/32/6 - 
- 2/22/42/6 



2/12/5 



+ 2/12/32/6 
2/12/42/6 



2/I2/6 



-2/4 + 22/32/42/5 

The coefficients ^4 through D are minors of a Hankel matrix and so Eq. (|A21|) can also be written compactly as 



2/1 


2/2 


2/3 


1 


2/2 


2/3 


2/4 




2/3 


2/4 


2/5 




2/4 


2/5 


2/6 


4 



The cubic equation can be solved using the method of Scipione del Ferro and Tartaglia [48 



General solution for M = 4 



(A22) 



Defining a fourth auxiliary quantity: 



the reduced system of equations for M 



4 is; 



Xiji^i --f^ {5<i< 2M) 



2/1 


— 01+02 


+ 03 + 04 


a2 


= 02 (xi - 


X2) + 03 (xi - X3) + 04 (xi - X4) 


Ps 


= 03 - 


X3) (X2 - .T3) + 04 {Xi - X4) (X2 - 


74 


= 04 (xi — 


X4,) {X2 - Xi) [X'i - Xi) 


^5 


= 




^6 


= 




S7 


= 




Ss 


= 





Xi) 



(A23) 



(A24) 
(A25) 
(A26) 
(A27) 
(A28) 
(A29) 
(A30) 
(A31) 



Following the now familiar procedure, substitute Eq. (|A23[) into the last four equations to eliminate x^ and find the 
(redundant) set of equations: 



7476 = 0, 7477 

7478 = 0, 7578 



- 7576 



0, ll 



- 1517 = 0, 



(A32) 



or as minors of a Hankel matrix: 



7677 = 0, 77 - 7678 = 



0, 



74 


75 


= 0, 


74 


75 


= 0, 


75 


76 


75 


76 


76 


77 




76 


77 


74 


76 


= 0, 


75 


76 


= 0, 


76 


77 


76 


78 


77 


78 


77 


78 



(A33) 



= 0. 



Substitute Eq. (|A11|) and eliminate X3 to find the next set of (redundant) set of equations: 

(3l - 2/34/35/36 + PsPl + PlPi - /33/35/37 



-/Jf/^e 
/35/3| 



PlPi + /33/35/38 



- 2/34/35/36 +/33/36 
■/34/35/37-/33/36/37 

• /34/36/37 + /33/3? + /34/35/38 - /33/36/38 

-al + 2/35/36/37 - /34/3? - /3|/38 + /34/36/38 



(A34) 
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or as minors of a Hankel matrix: 



/33 


(34 


(35 




(3z 


(3i 


(35 




(33 


(34 


(35 




04 


(35 


(36 






(3& 


= 0, 


(3i 


(35 


(36 




(35 


(36 


(37 


= 0, 


(35 


(36 


P7 


/35 




f3i 




(36 


f3i 


(3s 




(36 


(37 






(36 


(37 





= 0. 



(A35) 



Again, note that the LHS of these four equations are the same as the four coefficients of Eq. (|A21[) after the same 
substitution yi —i- Next, substitute Eq. (jA2p and eUminate X2 to find the equation 



2 2 2 2 2 

3 2 2 2 2 

'2a20i60i70i'^ + 2a3a4a8a5 — a2Q;g + a^a^ + agay 
-OL20L4oi^ — 2a3a4a60:7 — a^a^ — a^aeag + a20:4a6Cts = 



As in Eqs. (|A7ll and (|A20p the l.h.s. can be written as a determinant of a Hankel matrix of ai's: 



a2 as a4 as 

as a4 as ae 

a4 as ae ay 

as ae ay as 



FinaUy, substituting Eq. (jAip produces a quartic equation in xi 

= 



Vxl + Bx^J + (Jxl + i'xi + i!/' = U 

^ ^yt - iys.y5yl - 2y2y6yl - ymyl + '2y2yly4 + 2^32/62/4 + 2yiys2/62/' 
+ 2y2y3y7y4 - yiyf + + yhl - ymyl - 2^2^31/5^6 - yfy? 

- ^22/52/7 + 2/12/32/52/7 

- 2ys2/62/| - 2/4y7y| + 2y4yly3 + ^2^62/3 + 2/42/62/3 + 2/22/52/72/3 
+ 2/12/62/72/3 - 2y22/42/82/3 - 2/i2/52/82/3 - 2/22/5 " 2/i2/42/6 " 2/42/5 
+ 2/12/52/6 + 2/22/42/7 - 2/12/42/52/7 - 2/22/62/7 + 2/i2/42/8 + 2/22/52/8 

3 I 2 2 I 2 I 2 

=-2/62/4 + 2/52/4 + 2/32/72/4 

2 

■ 782/4 - 2/12/52/82/4 - !/3i/s ~ i/li/5i/6 ^2^7 ~ i/li/3i/7 "T 
'52/7 + 2/32/52/7 - 2/22/32/62/7 + 2/22/32/52/8 - 2/22/62/8 + 2/1 

+ 2/12/72/4 + 2/22/62/72/4 

' 2/32/52/6 + 2/22/^2/7 

+ 2/12/62/7 

22/32/62/5 ■ 

- 2/42/8 



-B =y82/| 



- 2/^2/82/4 - 2/12/52/82/4 - 2/32/f ^.^^^o ' ^^^y 

+ 2/12/52/7 + 2/32/52/7 - 2/22/32/62/7 + 2/22/32/52/8 - 

^ =-2/72/4 + 22/52/62/4 + 2/32/82/4 - 2/52/4 - 2?/32/62/4 + 2/12 

- 2/22/52/82/4 - 2/12/62/82/4 + 2/12/6 " 2/22/52/6 - 2/22/32/7 + 2/32/52/6 + 2/22/52/7 

+ 2/32/62/7 - 2?/iys2/62/7 + 2/i2/52/8 - 2/32/52/8 + 2/22/32/62/8 

^ =2/5 - 32/42/62/5 - 22/32/72/5 " 2/22/82/5 + 22/32/62/5 + 2^42/72/5 + 2J/22/62/72/5 

+ 2^32/42/82/5 - 2/22/6 + 2/42/6 + 2/32/7 - 2/22/42/7 - 2y32/42/62/7 

- 2/32/62/8 + 2/22/42/62/8 

As before, the coefficients A though i? are minors of a Hankel matrix of 2/„'s, so this equation can be written: 



2/22/82/4 + 2/22/62/4 - 32/22/52/72/4 

2,2 2 2 , 2 

' 2/22/52/6 



- 2/12/52/6 + 2/I2/7 - 2/12/32/7 

- 2/|y62/8 + 2/12/32/62/8 



(A36) 



(A37) 



(A38) 



2/1 


2/2 


2/3 


2/4 


1 


2/2 


2/3 


2/4 


2/5 




2/3 


2/4 


2/5 


2/6 


^? 


2/4 


2/5 


2/6 


2/7 




2/5 


2/6 


2/7 


2/8 


xt 



(A39) 



The quartic equation can be solved using the method of Ferrari [48 
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APPENDIX B: SOLUTION FOR EXAMPLES OF {K,M,N) 
1. Solution for {K, M, N) = (2, 2, 3) 
The nonlinear equations to solve has a block structure: 



2/11 




1 


1 










2/12 




Xl 










an 


2/13 




9 










ai2 


2/21 








1 


1 




a2i 


2/22 










a;2 




. a22 


- 2/23 . 








4 







(Bl) 



Here, the indices for y^n, and akm are in the ranges l<k<K,l<m<M and 1 < rt < A^. To reduce the 
system we extend Eq. to block form with K identical blocks Lm (x) on the diagonal. The reduced equations are 



2/fci = flfci 
Pk3 - 



+ ak2 
- X2) afc2 

(1 < fc < 2) 



(B2) 
(B3) 
(B4) 



where we have added an additional index k to the auxiliary quantities defined in Eqs. (jAl[) and (|A2[) . Substituting 
for Pk3 in Eqs. (jB4| and eliminating 2:2 gives the equation: 



ai2 OL22 

ai3 a23 



= 0, 



(B5) 



where we've written the equation as a minor of some matrix, following our experience in Appendix |^ yet whose 
structure is not yet clear. Substituting for gives a quadratic equation in x\ in determinant form: 



2/11 

2/12 
2/13 



2/21 
2/22 

2/23 



1 

Xl 



= 0. 



(B6) 



2. Solution for {K, M, N) = (2, 4, 6) 



The reduced system of equations for K = 2 correlation functions measured on iV = 6 equally spaced time slices to 
be solved to extract model parameters for M = 4 states is: 



2/fei = flfci + ak2 + afe3 + afe4 

ak2 = ak2 {xi - X2) + ak3 (xi - X3) + aki (xi - X4) 

/3fc3 = Ofes (Xi - X3) (X2 - X3) + flM (Xi - X4) (X2 - X4) 

7m = afe4 (a;i - X4) (x2 - X4) (X3 - X4) 
45 - 

^6 = 

(1 < fc < 2) 

Substituting for Skn and eliminating X4 gives a set of equations in 7fc„ : 



= 0, 



714 


724 


= 0, 


714 


725 


= 0, 


715 


724 


715 


725 




715 


726 




716 


725 


715 


725 


= 0, 


714 


715 


= 0, 


724 


725 


716 


726 




715 


716 




725 


726 



(B7) 
(B8) 
(B9) 
(BIO) 
(Bll) 
(B12) 



(B13) 
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Substituting for jkn and eliminating X3 gives a set of equations in (3kn 



/3l3 Pl4 P23 

Pl4 /3l5 P24 
Pl5 Pw P25 



0, 



/?13 /3l4 P24 

(3l4 Pl5 P25 
Pl5 P16 P26 



/3l3 P23 /324 
Pl4 1^24 P25 
Pl5 P25 P26 



Substituting for /3fe„ and eliminating X2 gives an equation in afe„: 



ai2 Q!13 a22 Q!23 

ai3 ai4 a23 024 

ai4 Q;i5 Q;24 "25 

Q15 aie a25 ce26 



Substituting for akn gives a quartic equation in Xi: 



yii yi2 


2/21 


2/22 


1 


yi2 2/13 


2/22 


2/23 




2/13 2/14 


2/23 


2/24 


xl 


2/14 2/15 


2/24 


2/25 


xl 


2/15 2/16 


2/25 


2/26 


t'^ 



f3l4 p23 P24 

/3l5 1^24 (^25 
P16 P25 P26 



(B14) 



(B15) 



(B16) 



3. Solution for (K,M,N) = (3,3,4) 



The reduced system of equations for iiT = 3 correlation functions measured on N = 4 equally spaced time slices to 
be solved to extract model parameters for M = 3 states is: 



2/fei — flfei + afe2 + dks 

otk2 = ak2 {xi - X2) + a/fcs {xi - xs) 

Pk3 = ak3 (Xl - X3) {X2 - X3) 
lk4 = 

(1 < fc < 3) 

Substituting for 7/54 and eliminating X3 gives a set of equations in 



013 


023 


= 0, 


013 


033 


= 0, 


023 


033 


014 


024 


014 


034 


024 


034 



Substituting for 0kn and eliminating X2 gives an equation in akn'- 



ai2 0.22 Ot32 
0L13 0.23 OL33 
ai4 a24 Q!34 



0. 



Substituting for akn gives a cubic equation in xi: 



yii 


2/21 


y31 


1 


2/12 


2/22 


2/32 


Xl 


2/13 


2/23 


2/33 


xl 


2/14 


2/24 


2/34 


xl 



0. 



0. 



(B17) 
(B18) 
(B19) 
(B20) 



(B21) 



(B22) 



(B23) 
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4. Solution for {K, M, N) = (4, 4, 5) 

The reduced system of equations ior K = A correlation functions measured on A/' = 5 equally spaced time slices to 
be solved to extract model parameters for M = 4 states is: 



Vki = aki + ak2 + cifes + ak4 (B24) 

Q.k2 = ak2 {xi - X2) + ak3 {xi - X3) + UkA {xi - X4) (B25) 

/3fc3 = ak3{xi-X3){x2-xa) + ak4{xi-X4){x2-X4) (B26) 

7/c4 = a/c4 (xi - X4) {x2 - X4) {x3 - X4) (B27) 

Sk5 = (B28) 
(1 < A; < 4) 

The determinant form of the quartic equation in Xi is 



yii 


2/21 


2/31 


2/41 


1 


yi2 


2/22 


2/32 


2/42 


Xl 


yi3 


2/23 


2/33 


2/43 


xl 


2/14 


2/24 


2/34 


2/44 


x\ 


2/15 


2/25 


2/35 


2/45 


x\ 
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